Frictional fluid instabilities shaped by viscous forces

Multiphase flows involving granular materials are complex and prone to pattern formation caused by competing mechanical and hydrodynamic interactions. Here we study the interplay between granular bulldozing and the stabilising effect of viscous pressure gradients in the invading fluid. Injection of aqueous solutions into layers of dry, hydrophobic grains represent a viscously stable scenario where we observe a transition from growth of a single frictional finger to simultaneous growth of multiple fingers as viscous forces are increased. The pattern is made more compact by the internal viscous pressure gradient, ultimately resulting in a fully stabilised front of frictional fingers advancing as a radial spoke pattern.

Multiphase flows involving granular materials are complex and prone to pattern formation caused by competing mechanical and hydrodynamic interactions. Here we study the interplay between granular bulldozing and the stabilising effect of viscous pressure gradients in the invading fluid. Injection of aqueous solutions into layers of dry, hydrophobic grains represent a viscously stable scenario where we observe a transition from growth of a single frictional finger to simultaneous growth of multiple fingers as viscous forces are increased. The pattern is made more compact by the internal viscous pressure gradient, ultimately resulting in a fully stabilised front of frictional fingers advancing as a radial spoke pattern.
Multiphase frictional flows inhabit a large parameter space, but relatively few scenarios have attracted any attention. One such flow that is now relatively well understood is the injection of a low-viscosity invading fluid (viscosity η inv ) to displace a much more viscous defending fluid (viscosity η def ≫ η inv ) containing sedimented grains, for the case where the invading fluid is nonwetting to the grains (i.e. drainage). Without the grains, or with grains that are fixed in place (i.e. within a rigid porous medium), this problem is famously viscously unstable and will be subject to classical viscous fingering (i.e. the Saffman-Taylor instability) [32][33][34][35] . With movable grains, the nonwetting invading phase will tend to bulldoze the defending mixture rather than invading the space between or above the grains as long as capillary forces are strong enough to overcome friction with the wall(s) and among the grains (i.e. the capillary entry pressure must be sufficiently larger than the frictional resistance to sliding and rearrangement). This bulldozing behaviour further destabilises the system as the accumulation of grains on the defending side of the interface penalises uniform displacement, leading to the formation of fractures, fingers, bubbles, labyrinths, and other patterns, depending on the injection rate and the packing fraction 15,19,28 .
Without the grains, reversing the two viscosities (η def < η inv ) negates the fingering instability by turning viscous pressure into a stabilising mechanism. With grains that are fixed in place, capillary forces and porescale disorder compete with viscous stabilisation to produce fractal invasion-percolation patterns at low injection rates and rough but stable fronts at high injection rates [34][35][36][37][38] . With movable grains, however, the flow is frictionally unstable at all rates due to bulldozing. The competition between these two mechanisms has not previously been studied in any detail and even basic questions remain unanswered; for example, to what extent can viscous forces stabilise the flow against the frictional instability? Here, we explore this competition systematically using experiments and simulations. Focusing on the case where η inv ≫ η def , we show that the pattern formation is controlled by the strength of viscous forces in the invading phase relative to friction due to bulldozing and pile-up of grains in the defending phase, as quantified by a "viscous deformability" parameter D visc . Increasing D visc leads to a striking transition from the growth of one solitary finger to the simultaneous growth of multiple, wandering fingers to the axisymmetric growth of a radial spoke pattern as the flow is increasingly viscously stabilised.

Viscously stable frictional fingers
Viscously stable frictional fingering was achieved by injecting water or a mixture of water and glycerol at a flow rate Q into a Hele-Shaw cell comprising two parallel glass plates separated by a gap of thickness b = 0.9 mm and containing a dry layer of polydisperse hydrophobic beads (Fig. 1a). As such, air was the low viscosity defending fluid in which the grains were "submerged". We define the log viscosity ratio M = logðη def =η inv Þ, which is negative for viscously stable scenarios and strongly negative here, where η def ≪ η inv . The flow cells were prepared with various filling levels φ = h/b of beads, where h is the initial thickness of the layer; the beads were silanised glass, sieved between 70 and 100 μm, with mean diameter d = 87 μm.
In all cases, the frictional instability shaped the invading water into one or more fingers of width 2R surrounded by a compaction front of thickness L of dry, bulldozed grains ( Fig. 1b and c). These fingers grew only at their tips; the side-walls were immobile after their initial formation, except where new fingers were initiated at larger values of Q, as discussed below. Figure 2a shows a selection of experimental results for different Q and φ. Figure 2b shows corresponding simulations described in a later section. Low injection rates resulted in the growth of a single finger that "wormed" its way through the cell (e.g. Q = 1 ml/min, bottom row in Fig. 2, Supplementary Movie 1). At the same injection rate, increasing φ led to a narrower finger because the compaction front thickened more rapidly with the advance of the interface (Fig. 3) 15,39 . Values for R were obtained from final images (such as those presented in Fig. 2) by measuring the ratio of fluid-filled invaded area to internal finger interface length R = A fluid /S finger .
The grains in the compaction front bridge the gap between the plates. In a straight segment of the side-wall, the capillary pressure imparted by the meniscus is opposed by the effective stress in the granular material, which disperses through grain-grain contacts to the plates. Deformation of the front requires dilation, which is opposed by the confining plates and the granular pressure from neighbouring front segments. Assuming Coulomb friction between the granular material and the plates and that out-of-plane stresses are proportional to the imposed in-plane stress (the "Janssen law" 40 ), the frictional stress resisting motion of the front increases exponentially with front width: σ(L) ∝ e L/ξ , where ξ = b/(2μκ) is the characteristic length, μ is the effective coefficient of friction between the grains and the plates, and κ is the ratio of out-of-plane normal stress to in-plane normal stress in the granular packing (i.e. the Janssen coefficient) 39,41 .
At the tip of a growing finger, the front is curved (Fig. 1c) and the streamlines of granular motion diverge as the grains are pushed outwards normal to the interface. In the reference frame of the moving tip, the granular material in the front is therefore continuously being stretched tangentially to the interface. This extensional motion reduces the granular pressure against the plates by weakening tangential granular force chains. The confinement-induced jamming that produced exponentially increasing friction at the straight side-walls is therefore significantly reduced at the moving tip, where the finger width is set. In agreement with previous work 42 , we find that a simple linear friction model provides a good fit to the experimental data in Fig. 3 (see model results described below). Following 42 , we estimate the frictional stress at the tip as σ t ≈ σ 0 L t /b, where L t is the front width at the tip (see Fig. 1c) and σ 0 represents the total frictional force per unit contact area with the plates; we use σ 0 as a fitting constant.
At the curved tip, the menisci in the pores between the grains collectively generate an effective interfacial tension at the scale of several grains that acts to minimise the curvature of the tip. The fluid pressure at the tip of a growing finger is therefore partly opposed by the Laplace pressure γ eff /R t associated with its in-plane curvature R t (Fig. 1), where γ eff is the effective interfacial tension at the liquid-gas/ grain interface 39 . We approximate the latter with the liquid-gas interfacial tension, γ eff ≈ γ. Note that we ignore the out-of-plane curvature because it is independent of the finger shape.
Following previous work 15,42,43 , we now estimate the characteristic rate-independent finger width 2R f that balances capillarity with friction by seeking the value of R f that minimises the capillary pressure at the tip P t . Conservation of mass for the finger as a whole suggests that (1 − φ)L f ≈ φR f , where L f is the characteristic front width and where we have neglected the small tip region; we assume for simplicity that the same argument applies independently at the tip, (1 − φ)L t ≈ φR t . From the arguments above, the total yield pressure at the fingertip is then suggesting that friction favours narrower fingers (displacing fewer grains) while surface tension favours wider fingers (less curvature). We then identify the R t at which P t is minimised by setting dP t /dR t = 0. Finally, we link R t to R f by requiring that this same yield pressure P t must apply along the straight side walls of the finger where the frictional stress is σ 0 L f /b and the curvature is negligible, noting that we are neglecting viscosity over the distance~2R f over which fingers transition from the curved tip to the straight side walls. Combining these ingredients, we arrive at a simple expression for R f : The white region has been invaded with water and cleared of grains, while the black region (the compaction front) has been completely filled with grains. c Advancing finger tips have radius of curvature R t and compaction-front width L t ; away from the tips, fully expanded fingers have half-width R and compaction-front width L. d Viscous flow of liquid along a growing finger leads to a pressure gradient from the injection pressure at the inlet to the capillary pressure P t at the tip.
This finger width 2R f is thus the emergent natural length scale in our system. A very similar expression has been derived and used for viscously unstable frictional fingers 15,42,43 .
We plot 2R as a function of φ in Fig. 3 for experiments and simulations at low rates Q, comparing against 2R f (φ) from Eq. (1). The theory agrees well with the experimental observations, in common with previous research on viscously unstable air-invasion labyrinths 15,42 and confirming the expectation that there exists a rateindependent regime where viscous forces within the fingers are indeed negligible, such that capillarity and friction compete to set the finger width.
Transition from a single finger to multiple fingers As we increase Q over two orders of magnitude, we observe a transition toward a regime where viscous pressure gradients become significant. As Q increases, the number of simultaneously growing fingers increases monotonically from one at low Q to nearly 20 at the highest values of Q explored here (Figs. 2 and 4, Supplementary Movie 2). Thus, strikingly, viscous forces drive the pattern to be more space-filling in this system by promoting the formation of more fingers.
Low-and high-Q experiments are compared in Fig. 4a and b, where the invading-fluid-filled fingers are colourised by invasion time. At low Q, growth is localised to the tip of a single finger into which all the injected fluid flows. Branching is sometimes observed, but usually only one finger is actively growing at any time. At high Q, several fingers grow simultaneously, as indicated by the rough axisymmetry of the colours in Fig. 4b. New fingers sprout by side-branching as the injected fluid flows through the network of fingers toward the active tips.
The viscous pressure gradient is located within the viscous invading fluid, with the pressure decreasing along the fingers from the inlet toward the moving tips (Fig. 1d). A sufficiently high pressure along the length of a finger can drive new fingers to break out from the side walls, leading to growth in the central parts of the pattern; this is the manifestation of viscous stabilisation in this frictionally unstable system. Figure 5a shows that the number of active fingers N increases with injected volume V and with injection rate Q, which is consistent with the fact that the fluid pressure at any fixed radial position along a finger must increase with the length of that finger and with the flow rate through it.
Breakout of new fingers is suppressed by at least two mechanisms. First, the stress required to deform a straight side-wall is higher than at the already curved tip, where the divergent flow of grains relieves bridging stresses by accommodating dilation. Second, the coefficient of static friction μ s at the side walls is likely to be higher than the coefficient of dynamic friction μ d at the moving tips, leading to higher frictional strength along the side walls. We represent these effects by introducing a higher threshold pressure P b = P t + ΔP b required to sprout a new finger from a side wall.
The results in §II A suggest that larger values of φ produce thicker compaction fronts, L f / ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi φ=ð1 À φÞ p 15 . These thicker fronts should provide a stronger resistance to breakout, meaning that ΔP b should increase with L f ; thus, it should become more difficult to form new   fingers as φ increases. This observation suggests that the number of fingers N should decrease with φ, as is confirmed in Fig. 5b. To capture this feature, we follow previous work 15,39,40 by hypothesising an exponential friction law of the form where we use the pre-factor σ β as a fitting parameter. The viscous contribution to the liquid pressure vanishes at a moving finger tip and increases linearly with distance upstream of the tip. The viscous pressure ΔP v a distance Δr b upstream of a moving finger tip can be estimated via Darcy's law, where Q/(2R f bN) is the average flux within each active finger. Thus, a new finger will form a distance Δr b behind a moving finger tip when ΔP v~Δ P b , suggesting that where r out is the radial system size (i.e. the outer radius) and is the dimensionless viscous deformability, which compares the characteristic viscous pressure drop to the characteristic frictional resistance of the side walls for a single finger of length r out . Note that Δr b decreases with D visc , so that faster injection or a more viscous invading phase will promote branching closer to the tip.
To sprout a single new finger from an existing one requires that a volume 2R f bΔr b of invading fluid be added to that finger. Thus, the addition of a volume ΔV to the flow cell will sprout ΔN new fingers, where For large N, we approximate the discrete variation ΔN as a continuous one and integrate from V = 0, giving where we have eliminated Δr b using Eq. (4). On the right-hand side of Eq. (7), the only unknown is ΔP b , which appears within D visc (recall that R f (φ) is given in Eq. (1)). We therefore plot the number of active fingers N against the total injected volume V from experiments at fixed φ for several different values of Q (Fig. 5a) and at fixed Q for several different values of φ (Fig. 5b). Then, we use ΔP b (φ) as a fitting parameter to achieve the best match between these results and the predictions of Eq. (7). In Fig. 5c, we then plot these bestfit values of ΔP b against ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi φ=ð1 À φÞ p / L f , as suggested by Eq. (2); the linear trend of the data on this semi-logarithmic scale provides support for the exponential friction model for the side-walls (Eq. (2)). A leastsquares fit of these empirical ΔP b values against ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi φ=ð1 À φÞ p (red line in Fig. 5c) suggests that σ β = 0.20 ± 0.02 Pa, where the 95% confidence interval of the fit is used for the uncertainty. Finally, we plot N(V) against for all of our experimental results in Fig. 5d (symbols), along with the prediction of Eq. (7) that these two quantities are directly proportional (red line).
The simple model captures the overall trend in the data, despite containing none of the geometrical complexity of the branching patterns. Note, however, that there is considerable variability in the data. The actual strength of viscous forces within a finger depends on the finger width, which exhibits 10-20% variability around the characteristic value R f at the same value of φ (Fig. 3). On the vertical axis, we observe that N can vary both between different experiments and also within a single experiment as fingers start, stop, and sometimes restart again, making this measurement inherently imprecise. Note also that, for large N, the fingers begin to fill the available space, crowding out the formation of new fingers. This effect is not included in the model, but should suppress the growth of N at even higher values of ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V D visc =ðR f br out Þ q . Nevertheless, our simple model captures reasonably well the roles of viscosity and friction in controlling the sprouting of new fingers.

Simulations of viscously stable frictional fingering
To better understand the complex patterns generated in experiments, we develop a frictional-fingering simulator by building on a code previously used to model air invasion into a wet hydrophilic packing 42,43 .
The granular filling-level field is represented as a 2D array (grey-scale pixels in Fig. 6a). The interface between the invading fluid and the defending mixture is represented by a chain of nodes (red points in the inset of Fig. 6a). As in the experiments, forward motion of the interface accumulates a compaction front. To implement viscous flow inside the invading phase, we skeletonise the finger structure into a branching tree (blue nodes and edges in the inset of Fig. 6), on which we calculate the local viscous pressure P v (Eq. (12)). Flow through the viscous skeleton is determined by the volume change dictated by growth at the tips of the fingers. An interface segment advances incrementally when its yield pressure is exceeded by the local fluid pressure, after which the filling fraction, the local interface curvatures R local , and the viscous pressures P v along the flow network are updated. Note that we neglect the viscosity of the defending fluid, and that we use two fixed fitting parameters, the friction coefficient and a viscous coefficient, both set to match finger widths and viscous stabilisation across the range of experimental parameters (see Methods). Figure 6a shows a detailed view of a simulation at moderate D visc = 22 with multiple active fingers. The skeletonised viscous flow network shows the fluid pressure, which decreases toward the tips of active fingers and is uniform along inactive fingers, with the maximum pressure at the inlet. Figure 6b shows the evolution of injection pressure for D visc = 31, comparing the experimental pressure measured at the inlet with the simulation pressure calculated at the central node for three realisations. In both experiment and simulations, the pressure initially increases rapidly as bulldozing mobilises friction and capillary forces, then more slowly due to build-up of viscous pressure as fingers grow in length. The simulations reproduce the trend in the pressure data from experiments, as well as the typical frequency and magnitude of the fluctuations; however, the overall magnitude of the injection pressure is somewhat higher in the simulations than in the experiment. Figure 2b shows simulation results across a range of filling fractions and injection rates, reproducing the experimental transition from single-finger to multi-finger growth as a function of Q, although somewhat under-predicting the number of fingers at high Q. Figure 3 includes the simulated finger width 2R at low Q, which decreases a bit more steeply with φ than what is observed in the experiments.
These differences between experiment and simulation may be due to the fact that the simulations use an exponential friction model (Eq. (11)) along the entire the compaction front, regardless of curvature. Our experimental data indicates that friction along the finger side walls is indeed well captured with an exponential model (Fig. 5), but that friction at the finger tips is better captured by a linear friction model (Fig. 3). We have used the exponential model in the simulations to prioritise viscous stabilisation, which is controlled by side-wall friction. However, tip friction controls the finger width and the resistance to finger propagation, so it is not surprising that the simulations exhibit a steeper variation of 2R with φ and a larger injection pressure than observed in experiments. Resolving this discrepancy would require the use of curvaturedependent friction along the compaction front, which may be the subject of future work.

From individual fingers to radial spoke patterns
To further increase the strength of viscous stabilisation relative to capillarity and friction, we increase the viscosity of the invading fluid. Figure 7 shows the time evolution of glycerol injection at Q = 10 mL/min, producing extreme viscous stabilisation and fingers that radiate outward in an axisymmetric spoke pattern. Here, the frictional instability produces fingers and viscous stabilisation forces them to grow radially, with an axisymmetric viscous pressure field. The tips remain equidistant from the inlet, creating a circular displacement front with embedded radial streaks of packed grains. As the pattern expands over time, the fingers increase in number by splitting to populate the growing circumference while maintaining a constant characteristic finger width (see Supplementary Movie 3). Figure 8a shows an experimental η inv -Q phase diagram in a loglog plot. The viscosity of the invading fluid η inv (i.e. the fraction of glycerol) increases from bottom to top and the injection rate Q increases from left to right. The top-right corner is empty because of the force limitation on the pump. Figure 8b shows a corresponding η inv -Q phase diagram for the simulations over a wider range of Q.
In these phase diagrams, similar patterns fall along diagonal lines corresponding to constant values of the product η inv Q, which is a measure of the strength of viscous forces. We quantify the balance of viscous forces, which drive the motion of the grains, to friction, which resists the motion of the grains, via the dimensionless viscous deformability D visc introduced above (Eq. (5)), which captures precisely this balance.
Our D visc is similar in spirit to the large-capillary-number limit of the 'fracturing number' of Holtzman et al. 21 , where the motion of a granular solid is resisted by friction under confining stress; to the 'viscous fracturing number' of Carrillo and Bourg 44 , where the motion of a porous viscoplastic solid is resisted by a yield stress; and to the 'fluidisation number' of Campbell et al. 25 , where the motion of a granular material was resisted by friction due to the weight of the grains. In the present context, friction is instead controlled by bulldozing, pile-up, and bridging.   8) and (9)). c Definitions of front length S front (red curve) and displaced area A dis (invaded area plus compaction front, blue region), with r max the reach of the most advanced finger. d Pattern compactness c = A dis /A reach , (e) front instability number s = S front /S finger , and (f) finger width 2R normalised by the characteristic rate-independent finger width 2R f (Eq. (1)) as functions of D visc . D * visc and D ** visc in blue and orange lines, and theoretical lower c 0 indicated in (d). In all panels, φ = 0.49, b = 0.9 mm, and the outer radius is r out = 13.4 cm.
Viscous stabilisation relates to the two macroscopic length scales: the system size, represented by the cell radius r out , and the emergent finger width 2R f . At low D visc , a single finger may grow to reach the perimeter of the cell without spawning a new finger if Δr b is greater than the system size r out . Taking Δr b = r out and N = 1, Eq. (4) suggests that the transition between a single finger and multiple fingers should occur around which is ultimately the scenario that motivated our definition of D visc . Note that D visc is proportional to system size, so a larger system requires a smaller viscosity or a lower injection rate to avoid branching. The next transition, from multiple fingers to radial spokes, occurs when the critical distance Δr b becomes smaller than the characteristic finger width 2R f . Branching will then be suppressed by the presence of neighbouring fingers as the pattern becomes space-filling. These patterns are viscously stabilised in the sense that a finger that advances ahead of its neighbours will decelerate due to its larger internal pressure drop, leading to a spoke pattern where all fingers extend roughly the same distance r from the inlet. The number of fingers N spoke (r) in this limit is easily estimated from conservation of mass, N spoke (r) = πr(1 − φ)/R f . Taking Δr b = 2R f and N = N spoke (r out ), Eq. (4) suggests that the transition to spokes should occur around ð9Þ We indicate these transitional values in the η-Q phase diagrams (diagonal blue and orange lines, respectively in Fig. 8a and b) for fixed φ = 0.49, b = 0.9 mm, and r out = 13.4 cm, for which D ** visc = 860. In both cases, we obtain a reasonable match to the visual characteristics of the patterns, transitioning from one finger to several around D ? visc and from multiple fingers to space-filling radial spokes around D ** visc , although both transitions are clearly gradual. Note that the D ? visc and D ** visc lines (blue and orange, respectively) appear steeper in Fig. 8b than in Fig. 8a because the range of Q is wider in the former.
For a quantitative analysis of the spatial characteristics of these patterns, we define the 'pattern compactness' c = A dis =ðπr 2 max Þ, where the displaced area A dis includes both fingers and compaction fronts (but not undisturbed material) and r max is the radial extent of the most advanced finger (Fig. 8c). For the measurements presented here, r max = r out . We plot c against D visc in Fig. 8d, indicating the relevant values of D * visc and D ** visc by the blue and orange lines, respectively. The theoretical lower limit c 0 = 2R f =ðπr out ð1 À φÞÞ ≈ 0:04 corresponds to a single straight finger growing from the inlet to the edge. The compactness c increases as viscous stabilisation creates multiple fingers and a more compact pattern, approaching 1 for radial spoke patterns.
We define a 'front instability number' s = S front /S finger , where S front is the length of the outer edge of the compaction front (i.e. the outer boundary between the compaction front and the undisturbed material; see red contour in Fig. 8c). For the spoke pattern, this boundary becomes nearly circular since the compaction fronts of neighbouring fingers touch. S finger is the longer internal perimeter of the finger pattern, tracing the liquid-air interface. The value of s is close to 1 for a single finger where the outer front perimeter follows the internal finger interface, decreasing as fingers increasingly meet to form a common front (Fig. 8e).
The finger width 2R (Fig. 8f) is approximately 2R f in the rateindependent regime (Eq. (1) and Fig. 3). Naively, one would expect increased viscous pressure within the finger to expand the width; instead, the fingers are observed to narrow when approaching D ** visc . This narrowing is most likely a result of self-confinement, in which the competition between numerous fingers increasingly suppresses lateral expansion. There is a gradual transition from multiple individual fingers to side-by-side radial spokes as viscous stabilisation becomes stronger and stronger, as evidenced by 2R beginning to decrease before the system reaches D ** visc . Note that this self-confinement effect is not included in the model.

Discussion
We have studied the fluid dynamics of an invading fluid displacing a defending fluid containing a sedimented granular material that is wetted by the defending fluid and repelled by the invading meniscus (drainage). For loose packings of grains, as in this study, capillary forces dominate over the weak frictional strength G of the initial granular layer, G = μ 0 ρ b gbφ, where μ 0 is the friction coefficient between the initial layer of grains and the plate, g is the body force per unit mass due to gravity, and ρ b = (ρ g − ρ def )(1 − n) is the bulk density difference between the granular layer and the defending fluid, where n is the porosity of the packing. In other words, the "capillary deformability" D cap = (γ/d)/G of the system is large in all experiments presented here, D cap ≫ 1, such that the meniscus can easily bulldoze the grains into compaction fronts that are then frictionally unstable, creating frictional fingers.
Holding capillarity constant, we then explored the effect of viscous stabilisation using mixtures of water and glycerol as the invading fluid and air as the defending fluid (negative M). The log viscosity contrast is large in all of our experiments (1:7 < |M| < 4:9), such that pressure gradients in the low-viscosity defending fluid (the air) are negligible.
For low Q and η inv , viscous forces are negligible on the scale of the finger width. Increasing the grain filling level φ increases the bulldozing friction and leads to narrower fingers. The finger width is set at the forward moving tip where a linear friction model σ(L f ) produces a good fit to the data. Breakout of new fingers from the static side-walls is suppressed by the frictional resistance of the granular compaction front, which increases exponentially as a function of its thickness.
Increasing Q or η inv increases the strength of viscous forces relative to frictional resistance, increasing the viscous deformability D visc . Viscously stable displacement involves pressure gradients along the invading fingers, with pressure decreasing from the central inlet toward the finger tips. Viscous stabilisation manifests as the sprouting of new fingers once the frictional 'breakout' pressure ΔP b of the walls is exceeded.
Two mechanisms determine the role of viscous stabilisation: (1) the strength of viscous pressure drop relative to frictional stress, as measured by D visc , and (2) the distance between the central inlet and the finger tips. We have identified two critical threshold values of D visc that separate different types of fingering patterns within the cell. Starting with a single finger at low D visc , increasing D visc eventually leads to the first threshold value at which the viscous pressure along the finger grows large enough to cause branching before the finger reaches the outer boundary; this value depends on the size of the system, since a longer finger implies a larger viscous pressure drop for the same value of D visc . As a result, larger cells would produce multiple fingers at lower D visc . Further increasing D visc leads to the second threshold value D ** visc , at which the viscous pressure gradient within the fingers is large enough to produce breakout pressures immediately behind the finger tips. Fingers that move ahead of the pack are suppressed by their internal viscous pressure drop and new fingers sprout continuously to populate an ever increasing pattern perimeter (in a radial cell). Ultimately, the fingers grow side-by-side in a space-filling radial spoke pattern. Figure 9 summarises the interplay between viscous stabilisation and friction in a D visc − φ phase diagram. Moving up the D visc axis increases number of active fingers and the compactness of the final pattern. Moving along the φ axis increases the frictional resistance to breakout, suppressing the viscous stabilisation mechanism. The estimated transitions from single to multiple fingers D * visc and from multiple fingers to spoke pattern D ** visc are plotted in blue and orange, respectively. Both transitions are gradual and the system exhibits a fairly large degree of variability, so neither value of D visc represents a sharp phase boundary; rather, they are indicative of the expected location of the transition region. As such, these transitional values agree reasonably well with the evolution of the geometric features ( Fig. 8d-f) and visual characteristics (Fig. 9) of the patterns.
Finally, it is instructive to consider our results in the context of traditional fluid displacement problems. As noted in the introduction, the corresponding fluid-fluid problem in a rigid Hele-Shaw cell is viscously stable and generates a circular displacement front in all cases. The defending phase in our system is a mixture of air and solid grains, so one might naively rationalise this frictional-fingering instability by thinking of the mixture of air and solid grains as a complex but effectively highly "viscous" defending phase, in which case the system would indeed be susceptible to a viscous-fingeringlike instability. However, this naive view is directly contradicted by the fact that our system is unstable at low rates and increasingly stabilised by viscosity at higher rates. In a rigid porous medium, this problem (viscously stable drainage) does have a pattern-forming mechanism at low rates that is stabilised by viscosity at higher rates, but those invasion-percolation patterns are fractal and controlled by pore-scale disorder; they are independent of any macroscopic physical properties and, by definition, lack a characteristic macroscopic length scale. In contrast, our problem features an emergent macroscopic finger-width that is much larger than the grain size and that varies in a predictable way with macroscopic interfacial tension, gap thickness, filling fraction, and frictional resistance (Eq. (1)). The fingering patterns themselves are weakly influenced by random spatial variations in the initial filling fraction, but insensitive to variations at the grain/pore scale.
Multiphase frictional flows are thus a distinct class of fluid displacement problems in which pattern formation is controlled by capillarity, viscosity, and both inter-granular and sliding friction. Drainage at large capillary deformability and strong mobility ratio is now relatively well understood for both stable and unstable scenarios (e.g. Ref. 19 and the present study, respectively), but much of the parameter space remains unexplored.

Experiments
The Hele-Shaw cell comprised two 40 × 40 × 1.5 cm glass plates separated by a gap thickness b = 0.9 mm. A 6 mm diameter hole through the centre of the top plate provided an inlet. The invading fluid was injected at controlled volume flow rates Q between 0.3 and 200 mL/ min using a syringe pump (Harvard Scientific, PHD Ultra). The cell was back-lit, and images were recorded using a Nikon 1 J2 digital camera at 30 fps.
The internal surfaces of the cell and the granular material were rendered hydrophobic by a silanization procedure following 45 . The silanization solution was a mixture of Trimethoxy(octadecyl)silane (OTMS) and Isopropyl alcohol(IPA). The silanization process was as follows: (1) the OTMS and IPA was mixed together in the ratio of 1:100; (2) the pH of the solution was adjusted to 3 by adding diluted Sulfuric acid (H 2 SO 4 , 0.1 M) to promote the hydrolysis of OTMS; (3) the solution was stirred using magnetic stirrer for at least 60 min at room temperature to form a alkylsilanol solution.
The glass pate surfaces were treated by pouring alkylsilanol solution on the surface and wiping over several times to make the coating uniform. The subsequently dried glass plates were hydrophobic with an estimated air/water contact angle of 120°. Silanization procedures were performed inside a fume hood.
The granular material was made hydrophobic by silanization treatment of soda-lime glass beads (Honite 18). The glass beads were acid-cleaned prior to silanization by the following steps: first, glass beads were immersed in hydrochloric acid (HCL, 0.1 M) and stirred using magnetic stirrer for at least 1 h. Then, they were rinsed thoroughly with deionized water and oven dried at 80°C. The dried beads were then sieved to a diameter range of 75-100 μm. The sieved beads were immersed into the silanization solution in a beaker and heated on a hotplate to accelerate the evaporation of the solution. The coated hydrophobic beads were sieved again to ensure no beads were clumped together.
Preparation of the granular layer: The dry hydrophobic beads were spread out on one of the treated glass plates (later to form the bottom surface of the Hele-Shaw cell). In order to achieve a granular layer of uniform thickness, two strips of adhesive tape were placed along opposite sides of the bottom plate, and a straight-edged tool resting on both tape strips was used to scrape the granular material into a uniform layer along the plate. The top plate was then mounted on top, separated from the bottom plate with 0.9 mm spacers. We varied the packing height by changing the thickness of the tape strips. Each strip consisted of several layers of tape film attached on top of one another. The tape film thickness was 63.5 μm, and between 6 and 12 layers were used to create strips producing granular layer heights h between 0.38 and 0.76 mm, corresponding to φ between 0.42 ± 0.01 and 0.84 ± 0.01, with the values verified and uncertainty estimated from trials where layers were made and then the mass measured independently. Experiments with lower filling fraction (φ = 0.21) are shown for formation of spoke patterns (Fig. 9), but note that filling fractions below 0.42 were not included in quantitative analysis of finger widths because of practical problems achieving uniform layer thickness for the thinnest layers. The cell was clamped together firmly after assembly to prevent the top plate from lifting. All four edges of the cell were left open to the atmosphere. Different viscosities were achieved by mixing glycerol and deionized water, with glycerol volume percentages of 0%, 58%, 84% and 100% corresponding to viscosities of 1, 14.14, 141.4 and 1414 mPa ⋅ s respectively 46 .

Simulations
The numerical code builds on the frictional fingering simulation presented in 39,43 . This code uses a two-dimensional array of values to represent the height of the sediment layer in the cell at each point in space, and a chain of nodes to represent the interface of the invading fluid, as illustrated in Fig. 6. Every timestep the modified threshold pressure P t is calculated for every interface node, following where R local is the local interface radius of curvature, σ is frictional stress resisting motion and P v is the viscous pressure difference between the inlet and that point. σ is expressed as in 15 : with L being the distance to the nearest point on the array which is not yet fully filled. The values used for bulk density was ρ b = 1450 kg/m 3 , and κ = 0.58 is the Janssen's coefficient. A static friction coefficient μ s = 0.91 is used, with the dynamic friction coefficient μ d = 0.9 being substituted in the exponent if the interface node has moved in the previous 500 cycles. The viscous pressure difference P v between the inlet and each point is calculated each timestep on a simplified version of the finger pattern, reduced to a branching tree of nodes as illustrated in Fig. 6. New nodes are added dynamically to this tree during the simulation whenever any section of interface is deemed to be too far away from its nearest node, to ensure that the shape of the tree closely mimics the shape of the invasion pattern. Each interface node reads P v from its nearest node on this tree. Each node calculates P v according to the Hagen-Poiseuille equation as where P 0 v is the P v of its parent node downstream, η is the viscosity, Q f is the flow rate into the finger downstream (averaged over the last 250 timesteps), X is the distance to its parent node and R is the average halfwidth of the finger between itself and its parent node. We use a calibration factor C = 1.33 to match the transitions in the simulations to the experimental observations. To calculate R, a two-dimensional array holds a value x for each position in the cell, x being the distance to the nearest interface; this array is updated whenever the interface advances. R is then estimated by stepping backwards towards the parent node along the ridge of the x distribution, taking the mean x along that ridge as the value of R.
Once the node with the lowest P t is identified, it is advanced forward slightly; its three nearest neighbours on each side are also moved by a lesser distance to maintain a smooth interface. New nodes are interpolated into the interface chain when the spacing between nodes exceeds a critical threshold, to maintain interface resolution as the interface lengthens. Whenever the interface advances, all granular material from the invaded region is redistributed to the nearest uninvaded positions which have space available. Randomness is introduced by initialising the distribution of granular material with random fluctuations above and below φ.
The 250 timestep period for averaging flow rate and the 500 timestep period for transitioning from kinetic to static friction are arbitrary numbers. They were chosen for being large compared to the typical number of growing fingers (to avoid excessive discretisation of viscous pressure and to prevent slower-moving but still active fingers from taking on static friction), but still small compared to the typical time period of an entire simulation. Qualitative tests did not suggest that the fingering patterns were strongly sensitive to these parameters. The simulated patterns match the experiments reasonably well across a wide range of D visc and φ.

Data availability
Images from experiments and simulations are available on the Zenodo data repository with the https://doi.org/10.5281/zenodo.7890690 47 .

Code availability
The Python code for the simulations are available from the authors upon request.